stanEx1
options(scipen=100,digits=3)
library(cmdstanr)
options(mc.cores=parallel::detectCores())
options(cmdstanr_max_rows=1000)
library(gridExtra)execute mcmc sampling
goMCMC=function(mdl,data,smp=500,wrm=100,th=1){
mcmc=mdl$sample(
data=data,
seed=1,
chains=4,
iter_sampling=smp,
iter_warmup=wrm,
thin=th,
refresh=1000
)
mcmc
}see mcmc result and parameters
seeMCMC=function(mcmc,exc='',ch=T){ # exclude 'exc' parameters from seeing
print(mcmc)
prs=mcmc$metadata()$model_params[-1] # reject lp__
for(pr in prs){
if(grepl('^y',pr)) next # not show predictive value "y*" information
if(exc!='' && grepl(paste0('^',exc),pr)) next
drw=mcmc$draws(pr)
if(ch){
par(mfrow=c(4,1),mar=c(1,5,1,1))
drw[,1,] |> plot(type='l',xlab='',ylab=pr)
drw[,2,] |> plot(type='l',xlab='',ylab=pr)
drw[,3,] |> plot(type='l',xlab='',ylab=pr)
drw[,4,] |> plot(type='l',xlab='',ylab=pr)
par(mar=c(3,5,3,3))
}
par(mfrow=c(1,2))
drw |> hist(main=pr,xlab='')
drw |> density() |> plot(main=pr)
}
}stan basic
normal distribution
ex1-0.stan
normal distribution
data{
int N;
vector[N] y;
}
parameters{
real m;
real<lower=0> s;
}
model{
y~normal(m,s);
}
generated quantities{
real y1;
y1=normal_rng(m,s);
}mdl=cmdstan_model('./ex1-0.stan')
y=rnorm(10,2,1)
summary(y)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.16 1.48 2.25 2.21 2.93 4.20par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y') 
data=list(N=length(y),y=y)
mle=mdl$optimize(data=data)## Initial log joint probability = -25.8051
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 12 -6.3269 0.000301942 0.000252394 0.9807 0.9807 14
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.2 seconds.mle## variable estimate
## lp__ -6.33
## m 2.21
## s 1.14
## y1 2.91mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.4 seconds.mcmc## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -7.31 -6.97 1.13 0.83 -9.41 -6.23 1.00 593 596
## m 2.23 2.22 0.48 0.43 1.48 3.01 1.00 627 642
## s 1.43 1.34 0.42 0.35 0.93 2.22 1.00 1127 990
## y1 2.24 2.23 1.57 1.40 -0.30 4.77 1.00 1676 1643mcmc$metadata()$stan_variables## [1] "lp__" "m" "s" "y1"mcmc$metadata()$model_params## [1] "lp__" "m" "s" "y1"seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -7.31 -6.97 1.13 0.83 -9.41 -6.23 1.00 593 596
## m 2.23 2.22 0.48 0.43 1.48 3.01 1.00 627 642
## s 1.43 1.34 0.42 0.35 0.93 2.22 1.00 1127 990
## y1 2.24 2.23 1.57 1.40 -0.30 4.77 1.00 1676 1643
poisson distribution
y=rpois(20,3)
summary(y)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.00 2.00 2.50 2.85 4.00 6.00par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')
data=list(N=length(y),y=y)ex2-1.stan
in case fit poisson dist. to normal dist.
data{
int N;
vector[N] y;
}
parameters{
real<lower=0> m;
real<lower=0> s;
}
model{
y~normal(m,s);
}
generated quantities{
real y1;
y1=normal_rng(m,s);
}# when fitting poisson dist. sample to normal dist.
mdl=cmdstan_model('./ex2-1.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -174.991
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 20 -16.5622 2.48557e-05 0.000216293 1 1 26
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -16.56
## m 2.85
## s 1.39
## y1 1.96mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -16.22 -15.91 1.05 0.75 -18.41 -15.22 1.00 882 990
## m 2.84 2.84 0.34 0.33 2.28 3.40 1.00 1580 1197
## s 1.52 1.48 0.27 0.24 1.15 2.04 1.00 851 860
## y1 2.80 2.80 1.57 1.57 0.21 5.29 1.00 2134 2016
y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1') 
ex2-2.stan
fit to poisson dist.
data{
int N;
array[N] int y;
}
parameters{
real<lower=0> l;
}
model{
y~poisson(l);
}
generated quantities{
int y1;
y1=poisson_rng(l);
}mdl=cmdstan_model('./ex2-2.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -18.63
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 6 2.69718 0.00015116 4.20026e-05 1 1 8
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ 2.70
## l 2.85
## y1 0.00mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ 3.25 3.52 0.71 0.32 1.88 3.75 1.00 1039 1194
## l 2.91 2.90 0.38 0.38 2.28 3.56 1.00 735 1050
## y1 2.84 3.00 1.71 1.48 0.00 6.00 1.00 1833 1969
y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1') 
difference test
ex3.stan
mean difference of 2 normal distributions
data{
int N;
vector[N] a;
vector[N] b;
}
parameters{
real ma;
real<lower=0> sa;
real mb;
real<lower=0> sb;
}
model{
a~normal(ma,sa);
b~normal(mb,sb);
}
generated quantities{
real d;
d=mb-ma;
}n=20
a=rnorm(n,10,1)
b=rnorm(n,11,2)
data=list(N=n,a=a,b=b)
mdl=cmdstan_model('./ex3.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -13861
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 21 -35.6979 0.000837279 0.000622165 1 1 48
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -35.70
## ma 9.78
## sa 1.12
## mb 10.58
## sb 1.96
## d 0.80mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -36.98 -36.65 1.48 1.30 -39.68 -35.25 1.00 802 970
## ma 9.78 9.79 0.27 0.27 9.35 10.21 1.00 1611 1325
## sa 1.22 1.20 0.22 0.21 0.93 1.62 1.00 2402 1684
## mb 10.58 10.59 0.49 0.48 9.78 11.36 1.00 616 489
## sb 2.16 2.10 0.40 0.35 1.66 2.90 1.00 1729 1391
## d 0.80 0.80 0.56 0.53 -0.15 1.72 1.00 698 814
d=mcmc$draws('d')
d## # A draws_array: 500 iterations, 4 chains, and 1 variables
## , , variable = d
##
## chain
## iteration 1 2 3 4
## 1 0.66 0.53 0.42 0.49
## 2 -0.25 0.63 1.33 0.93
## 3 -0.36 0.60 1.14 0.83
## 4 0.68 0.66 0.90 0.53
## 5 0.45 0.67 1.32 1.07
##
## # ... with 495 more iterationsmean(d>0)## [1] 0.916normal regression
ex4-1.stan
single normal regression
data{
int N;
vector[N] x;
vector[N] y;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~normal(b0+b1*x,s);
}# make prediction with R
n=20
x=runif(n,0,100)
y=rnorm(n,x*3+10,10)
par(mfrow=c(1,1))
plot(x,y)
data=list(N=n,x=x,y=y)
mdl=cmdstan_model('./ex4-1.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -4.34345e+06
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/RtmpSt0wv7/model-a38248be373.stan', line 14, column 2 to column 22)
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## 29 -50.2516 0.000701789 0.000697943 1 1 84
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -50.25
## b0 10.52
## b1 3.04
## s 7.48mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -49.89 -49.49 1.45 1.13 -52.79 -48.42 1.01 310 306
## b0 10.32 10.38 3.86 3.41 3.71 16.62 1.04 105 85
## b1 3.04 3.04 0.08 0.07 2.92 3.16 1.03 139 132
## s 8.50 8.26 1.65 1.45 6.27 11.53 1.00 714 459
b0=mcmc$draws('b0')
b1=mcmc$draws('b1')
b=tibble(b0=c(b0),b1=c(b1)) |> sample_n(100)
x1=runif(nrow(b),min(x),max(x))
xy=tibble(x=x1,y=b$b0+b$b1*x1)
par(mfrow=c(1,1))
plot(xy)
ex4-2.stan
single normal regression, prediction from new explanatory variable x1
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x1;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~normal(b0+b1*x,s);
}
generated quantities{
vector[N1] m;
vector[N1] y1;
for(i in 1:N1){
m[i]=b0+b1*x1[i];
y1[i]=normal_rng(m[i],s);
}
}# make prediction with stan
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex4-2.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -1.89582e+06
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## 49 -50.2516 0.000264535 0.000143774 0.1514 1 104
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -50.25
## b0 10.52
## b1 3.04
## s 7.48
## m[1] 15.50
## m[2] 46.15
## m[3] 76.79
## m[4] 107.44
## m[5] 138.09
## m[6] 168.73
## m[7] 199.38
## m[8] 230.03
## m[9] 260.67
## m[10] 291.32
## y1[1] 13.75
## y1[2] 57.26
## y1[3] 72.47
## y1[4] 94.83
## y1[5] 132.97
## y1[6] 172.74
## y1[7] 198.18
## y1[8] 207.57
## y1[9] 257.93
## y1[10] 297.49mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.mcmc$metadata()$stan_variables## [1] "lp__" "b0" "b1" "s" "m" "y1"seeMCMC(mcmc,'m',ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -49.85 -49.49 1.33 1.05 -52.42 -48.42 1.01 322 298
## b0 10.59 10.51 4.03 3.69 4.37 17.69 1.02 214 143
## b1 3.04 3.04 0.08 0.07 2.91 3.16 1.01 266 224
## s 8.48 8.30 1.49 1.41 6.42 11.21 1.00 637 1083
## m[1] 15.57 15.49 3.92 3.61 9.50 22.48 1.02 214 143
## m[2] 46.20 46.20 3.27 3.04 41.10 51.87 1.02 216 134
## m[3] 76.83 76.87 2.69 2.54 72.64 81.46 1.02 228 185
## m[4] 107.46 107.48 2.23 2.14 104.00 111.04 1.01 277 340
## m[5] 138.10 138.09 1.97 1.93 134.89 141.37 1.01 491 881
## m[6] 168.73 168.71 2.00 1.97 165.42 172.01 1.00 1573 1313
## m[7] 199.36 199.31 2.30 2.34 195.56 203.15 1.00 2253 1537
## m[8] 230.00 230.02 2.79 2.76 225.35 234.58 1.00 1479 1243
## m[9] 260.63 260.60 3.39 3.32 255.17 266.10 1.00 885 740
## m[10] 291.26 291.27 4.05 3.96 284.87 297.78 1.00 615 604
## y1[1] 15.74 15.77 9.56 9.21 0.46 31.26 1.00 740 1342
## y1[2] 46.14 45.89 9.49 9.02 30.69 61.30 1.01 944 1562
## y1[3] 76.64 76.66 9.29 8.72 61.64 92.04 1.00 1463 1702
## y1[4] 107.67 107.82 8.80 8.62 93.34 121.44 1.00 1608 1785
## y1[5] 137.99 138.15 8.85 8.20 123.56 151.94 1.00 1622 1737
## y1[6] 168.59 168.67 8.92 8.58 153.36 183.19 1.00 1951 1892
## y1[7] 198.92 198.80 8.85 8.94 184.93 213.30 1.00 1938 1937
## y1[8] 230.17 230.19 9.34 9.02 214.18 245.01 1.00 2132 1964
## y1[9] 260.66 260.79 9.15 8.58 245.69 275.14 1.00 1913 1871
## y1[10] 291.48 291.27 9.44 9.16 276.06 307.20 1.00 1525 1648
m=mcmc$draws('m')
summary(m)## # A tibble: 10 × 10
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 m[1] 15.6 15.5 3.92 3.61 9.50 22.5 1.02 214. 143.
## 2 m[2] 46.2 46.2 3.27 3.04 41.1 51.9 1.02 216. 135.
## 3 m[3] 76.8 76.9 2.69 2.54 72.6 81.5 1.02 229. 186.
## 4 m[4] 107. 107. 2.23 2.14 104. 111. 1.01 278. 340.
## 5 m[5] 138. 138. 1.97 1.93 135. 141. 1.01 492. 882.
## 6 m[6] 169. 169. 2.00 1.97 165. 172. 1.00 1574. 1313.
## 7 m[7] 199. 199. 2.30 2.34 196. 203. 1.00 2253. 1538.
## 8 m[8] 230. 230. 2.79 2.76 225. 235. 1.00 1480. 1243.
## 9 m[9] 261. 261. 3.39 3.32 255. 266. 1.00 886. 741.
## 10 m[10] 291. 291. 4.05 3.96 285. 298. 1.00 615. 604.smm=summary(m)
y1=mcmc$draws('y1')
summary(y1)## # A tibble: 10 × 10
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 y1[1] 15.7 15.8 9.56 9.21 0.459 31.3 1.00 740. 1342.
## 2 y1[2] 46.1 45.9 9.49 9.02 30.7 61.3 1.01 945. 1563.
## 3 y1[3] 76.6 76.7 9.29 8.72 61.6 92.0 1.00 1463. 1702.
## 4 y1[4] 108. 108. 8.80 8.62 93.3 121. 1.00 1609. 1785.
## 5 y1[5] 138. 138. 8.85 8.20 124. 152. 1.00 1622. 1737.
## 6 y1[6] 169. 169. 8.92 8.58 153. 183. 1.00 1952. 1893.
## 7 y1[7] 199. 199. 8.85 8.94 185. 213. 0.999 1939. 1938.
## 8 y1[8] 230. 230. 9.34 9.02 214. 245. 1.00 2133. 1964.
## 9 y1[9] 261. 261. 9.15 8.58 246. 275. 1.00 1914. 1871.
## 10 y1[10] 291. 291. 9.44 9.16 276. 307. 1.00 1525. 1649.smy=summary(y1)
xy=tibble(x=x1,m=smm$median,ml5=smm$q5,mu5=smm$q95,yl5=smy$q5,yu5=smy$q95)
par(mfrow=c(1,1))
xlim=c(min(x),max(x))
ylim=c(min(y),max(y))
plot(x,y,
xlim=xlim,ylim=ylim, xlab='x',ylab='y',col='red')
par(new=T)
plot(xy$x,xy$m,type='l',
xlim=xlim,ylim=ylim, xlab='',ylab='',col='black')
par(new=T)
plot(xy$x,xy$ml5,type='l',
xlim=xlim,ylim=ylim, xlab='',ylab='',col='darkgray')
par(new=T)
plot(xy$x,xy$mu5,type='l',
xlim=xlim,ylim=ylim, xlab='',ylab='',col='darkgray')
par(new=T)
plot(xy$x,xy$yl5,type='l',
xlim=xlim,ylim=ylim, xlab='',ylab='',col='lightgray')
par(new=T)
plot(xy$x,xy$yu5,type='l',
xlim=xlim,ylim=ylim, xlab='',ylab='',col='lightgray')
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=m),xy,col='black')+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')
ex4-3.stan
single normal regression, prediction from original explanatory variable x
data{
int N;
vector[N] x;
vector[N] y;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~normal(b0+b1*x,s);
}
generated quantities{
vector[N] m;
vector[N] y1;
for(i in 1:N){
m[i]=b0+b1*x[i];
y1[i]=normal_rng(m[i],s);
}
}data=list(N=n,x=x,y=y)
mdl=cmdstan_model('./ex4-3.stan')
mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,'m',ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -49.84 -49.51 1.29 1.08 -52.42 -48.42 1.00 495 534
## b0 10.54 10.66 3.68 3.41 4.63 16.37 1.01 222 223
## b1 3.04 3.03 0.07 0.07 2.92 3.15 1.01 297 428
## s 8.47 8.29 1.51 1.37 6.39 11.30 1.00 1030 848
## m[1] 140.85 140.85 1.99 1.90 137.49 143.98 1.00 685 921
## m[2] 255.85 255.84 3.35 3.28 250.01 261.11 1.00 1072 1166
## m[3] 16.72 16.81 3.56 3.30 10.97 22.33 1.01 222 225
## m[4] 56.66 56.75 2.82 2.64 52.11 61.17 1.01 233 238
## m[5] 218.20 218.24 2.69 2.67 213.64 222.41 1.00 2060 1440
## m[6] 15.51 15.61 3.58 3.31 9.75 21.15 1.01 222 225
## m[7] 213.08 213.15 2.61 2.55 208.65 217.15 1.00 2168 1519
## m[8] 148.88 148.86 1.99 1.90 145.52 152.00 1.00 896 1146
## m[9] 83.91 84.02 2.41 2.29 79.98 87.80 1.01 264 313
## m[10] 97.86 97.92 2.24 2.13 94.18 101.39 1.01 296 433
## m[11] 236.23 236.24 2.99 2.97 230.98 240.92 1.00 1460 1311
## m[12] 72.22 72.34 2.58 2.43 67.98 76.37 1.01 247 256
## m[13] 93.37 93.43 2.29 2.17 89.62 97.01 1.01 284 385
## m[14] 88.86 88.93 2.34 2.23 85.04 92.59 1.01 273 332
## m[15] 103.85 103.91 2.18 2.10 100.26 107.25 1.01 316 443
## m[16] 210.47 210.53 2.57 2.50 206.12 214.48 1.00 2222 1520
## m[17] 291.17 291.18 4.05 3.95 284.29 297.59 1.00 754 949
## m[18] 203.46 203.50 2.47 2.40 199.30 207.34 1.00 2304 1496
## m[19] 130.83 130.81 2.00 1.89 127.45 133.94 1.00 518 756
## m[20] 229.93 229.96 2.88 2.86 224.98 234.45 1.00 1646 1370
## y1[1] 140.87 140.80 8.59 8.21 127.16 155.12 1.00 2031 1650
## y1[2] 255.86 255.66 9.48 9.19 240.35 271.57 1.00 1683 1732
## y1[3] 16.66 16.90 9.24 9.05 1.46 31.96 1.00 1080 1446
## y1[4] 56.75 56.84 9.22 8.88 42.54 72.00 1.00 1223 1788
## y1[5] 218.14 218.27 9.09 8.69 203.36 232.98 1.00 2180 1884
## y1[6] 15.22 15.52 9.31 8.86 -0.80 30.36 1.00 876 1357
## y1[7] 213.51 213.57 9.04 8.85 198.64 228.31 1.00 2142 1758
## y1[8] 148.90 148.88 8.92 8.48 133.98 163.38 1.00 1877 1869
## y1[9] 84.08 83.81 8.96 8.62 69.58 98.84 1.00 1555 1711
## y1[10] 97.75 97.65 8.90 8.85 83.16 112.47 1.00 1476 1645
## y1[11] 236.10 236.07 9.09 8.86 221.04 250.92 1.00 1678 1493
## y1[12] 71.99 72.23 9.05 8.63 57.11 86.79 1.00 1311 1522
## y1[13] 93.28 93.32 9.05 8.48 78.02 107.96 1.00 1558 1622
## y1[14] 88.88 88.96 8.87 8.33 74.14 103.85 1.00 1691 1813
## y1[15] 103.96 103.72 8.85 8.45 89.81 118.65 1.00 1886 1935
## y1[16] 210.08 210.37 8.79 8.49 195.21 224.62 1.00 1831 1749
## y1[17] 290.90 291.08 9.57 9.11 275.56 306.53 1.00 1760 1735
## y1[18] 203.51 203.47 9.01 8.67 188.48 217.93 1.00 1626 1892
## y1[19] 131.12 130.96 8.64 8.12 116.99 145.29 1.00 1539 1912
## y1[20] 229.91 230.18 8.99 8.41 215.21 244.50 1.00 1986 1892
y1=mcmc$draws('y1')
smy=summary(y1)
par(mfrow=c(1,1))
plot(y,smy$median,xlab='obs.',ylab='prd.')
abline(0,1)
qplot(y,smy$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1)
par(mfrow=c(1,2))
hist(y-smy$median,xlab='obs.-prd.',main='residual')
density(y-smy$median) |> plot(xlab='obs.-prd.',main='residual')
estimate correlation
ex4-41.stan
use covariance matrix and multinormal distribution
data{
int N;
array[N] vector[2] y;
}
parameters{
vector[2] m;
vector<lower=0>[2] s;
real<lower=-1,upper=1> r;
}
transformed parameters{
cov_matrix[2] cv;
cv[1,1]=s[1]^2;
cv[2,2]=s[2]^2;
cv[1,2]=r*s[1]*s[2];
cv[2,1]=r*s[1]*s[2];
}
model{
y~multi_normal(m,cv);
}ex4-42.stan
use covariance matrix and wishart distribution
data{
int N;
matrix[2,2] dp;
}
parameters{
vector<lower=0>[2] s;
real<lower=-1,upper=1> r;
}
transformed parameters{
cov_matrix[2] cv;
cv[1,1]=s[1]^2;
cv[2,2]=s[2]^2;
cv[1,2]=r*s[1]*s[2];
cv[2,1]=r*s[1]*s[2];
}
model{
dp~wishart(N-1,cv);
}ex4-43.stan
use correlation matrix and wishart distribution
data{
int N;
int M;
matrix[M,M] dp;
}
parameters{
vector<lower=0>[M] s;
corr_matrix[M] r;
}
transformed parameters{
cov_matrix[M] cv;
cv=quad_form_diag(r,s);
}
model{
dp~wishart(N-1,cv);
}ex4-44.stan use covariance matrix and multinormal distribution
data{
int N;
int K;
array[N] vector[K] y;
}
parameters{
vector[K] m;
cov_matrix[K] cov;
}
model{
y~multi_normal(m,cov);
}library(MASS)
n=20
y=mvrnorm(n,c(0,0),matrix(c(1,0.7,0.7,1),2),2)
cov(y)## [,1] [,2]
## [1,] 0.857 0.664
## [2,] 0.664 1.158cor(y)## [,1] [,2]
## [1,] 1.000 0.666
## [2,] 0.666 1.000plot(y)
#estimate covariance
mdl=cmdstan_model('./ex4-41.stan')
data=list(N=n,y=y)
mle=mdl$optimize(data=data)## Initial log joint probability = -120.143
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 22 -13.0271 3.32791e-05 5.95564e-05 1 1 23
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -13.03
## m[1] -0.18
## m[2] -0.18
## s[1] 0.90
## s[2] 1.05
## r 0.67
## cv[1,1] 0.81
## cv[2,1] 0.63
## cv[1,2] 0.63
## cv[2,2] 1.10mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -17.13 -16.73 1.81 1.64 -20.35 -14.95 1.01 699 1006
## m[1] -0.18 -0.19 0.24 0.23 -0.57 0.21 1.00 1137 1125
## m[2] -0.19 -0.20 0.28 0.26 -0.62 0.27 1.00 1101 968
## s[1] 1.01 0.98 0.18 0.17 0.75 1.34 1.00 1372 1342
## s[2] 1.17 1.13 0.21 0.19 0.88 1.56 1.00 1311 1465
## r 0.62 0.64 0.15 0.14 0.34 0.82 1.00 614 831
## cv[1,1] 1.05 0.96 0.41 0.32 0.56 1.79 1.00 1372 1342
## cv[2,1] 0.76 0.69 0.37 0.30 0.30 1.44 1.00 711 1031
## cv[1,2] 0.76 0.69 0.37 0.30 0.30 1.44 1.00 711 1031
## cv[2,2] 1.41 1.29 0.55 0.42 0.78 2.43 1.00 1311 1465
#estimate correlation,covariance
#deviation product sum matirix folows Wishart dist.
#deviation product sum = covariance * n-1 or n
mdl=cmdstan_model('./ex4-42.stan')
data=list(N=n,dp=cov(y)*(n-1))
mle=mdl$optimize(data=data)## Initial log joint probability = -52.4202
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 15 -13.3503 0.000292043 3.38327e-05 1 1 17
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -13.35
## s[1] 0.93
## s[2] 1.08
## r 0.67
## cv[1,1] 0.86
## cv[2,1] 0.66
## cv[1,2] 0.66
## cv[2,2] 1.16mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -16.17 -15.81 1.32 1.01 -18.79 -14.79 1.00 874 978
## s[1] 1.01 0.98 0.19 0.17 0.75 1.33 1.00 1088 1036
## s[2] 1.17 1.14 0.22 0.19 0.88 1.55 1.00 1114 1260
## r 0.62 0.64 0.14 0.14 0.37 0.82 1.00 500 693
## cv[1,1] 1.05 0.96 0.41 0.32 0.56 1.77 1.00 1088 1036
## cv[2,1] 0.76 0.70 0.38 0.30 0.33 1.49 1.00 581 793
## cv[1,2] 0.76 0.70 0.38 0.30 0.33 1.49 1.00 581 793
## cv[2,2] 1.41 1.29 0.57 0.42 0.78 2.41 1.00 1114 1260
#estimate correlation,covariance
mdl=cmdstan_model('./ex4-43.stan')
m=2
data=list(N=n,M=m,dp=cov(y)*(n-1))
mle=mdl$optimize(data=data)## Initial log joint probability = -36.1841
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 11 -13.3503 0.000599028 0.000839941 1 1 15
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -13.35
## s[1] 0.93
## s[2] 1.08
## r[1,1] 1.00
## r[2,1] 0.67
## r[1,2] 0.67
## r[2,2] 1.00
## cv[1,1] 0.86
## cv[2,1] 0.66
## cv[1,2] 0.66
## cv[2,2] 1.16mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -15.57 -15.23 1.36 1.10 -18.32 -14.10 1.00 775 1110
## s[1] 1.01 0.98 0.19 0.18 0.76 1.33 1.00 1220 906
## s[2] 1.17 1.15 0.22 0.20 0.88 1.57 1.00 946 1030
## r[1,1] 1.00 1.00 0.00 0.00 1.00 1.00 NA NA NA
## r[2,1] 0.61 0.64 0.15 0.14 0.32 0.83 1.00 741 891
## r[1,2] 0.61 0.64 0.15 0.14 0.32 0.83 1.00 741 891
## r[2,2] 1.00 1.00 0.00 0.00 1.00 1.00 NA NA NA
## cv[1,1] 1.05 0.97 0.43 0.34 0.57 1.78 1.00 1220 906
## cv[2,1] 0.76 0.68 0.41 0.32 0.29 1.52 1.00 719 663
## cv[1,2] 0.76 0.68 0.41 0.32 0.29 1.52 1.00 719 663
## cv[2,2] 1.43 1.32 0.58 0.45 0.77 2.47 1.00 946 1030
mdl=cmdstan_model('./ex4-44.stan')
k=2
data=list(N=n,K=k,y=y)
mle=mdl$optimize(data=data)## Initial log joint probability = -624.153
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 14 -13.0271 0.000339274 0.000683859 1 1 16
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -13.03
## m[1] -0.18
## m[2] -0.18
## cov[1,1] 0.81
## cov[2,1] 0.63
## cov[1,2] 0.63
## cov[2,2] 1.10mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -15.00 -14.65 1.74 1.58 -18.35 -12.85 1.00 761 933
## m[1] -0.18 -0.18 0.25 0.24 -0.60 0.23 1.00 1338 1237
## m[2] -0.18 -0.18 0.29 0.28 -0.64 0.28 1.00 1162 1339
## cov[1,1] 1.25 1.12 0.56 0.40 0.66 2.25 1.00 1444 1098
## cov[2,1] 0.97 0.86 0.53 0.37 0.35 1.93 1.00 1085 817
## cov[1,2] 0.97 0.86 0.53 0.37 0.35 1.93 1.00 1085 817
## cov[2,2] 1.67 1.50 0.70 0.49 0.90 3.00 1.00 1223 850
distributions
normal dist.
distribution of mean, it's from central limit theorem
y~N(m,s), y(-Inf,Inf), E[y]=m, V[y]=s^2ex1-0.stan
data{
int N;
vector[N] y;
}
parameters{
real m;
real<lower=0> s;
}
model{
y~normal(m,s);
}
generated quantities{
real y1;
y1=normal_rng(m,s);
}n=20
y=rnorm(n,2,1)
summary(y)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.53 1.45 1.98 1.98 2.47 3.21par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')
data=list(N=n,y=y)
mdl=cmdstan_model('./ex1-0.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -60.248
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 13 -3.5906 3.72627e-05 4.97546e-05 1 1 14
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -3.59
## m 1.98
## s 0.73
## y1 2.96mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -4.95 -4.60 1.11 0.70 -7.22 -3.95 1.00 762 920
## m 1.98 1.99 0.18 0.17 1.70 2.26 1.00 1628 1082
## s 0.80 0.78 0.14 0.13 0.61 1.07 1.00 1682 1166
## y1 1.96 1.97 0.86 0.84 0.58 3.34 1.00 1999 1932
Bernoulli dist.
The event with probablity p occur y=1, not occur y=0
y~Ber(p), y{0,1}, E[y]=p, V[y]=p(1-p)ex1-1.stan
data{
int N;
array[N] int y;
}
parameters{
real<lower=0,upper=1> p;
}
model{
y~bernoulli(p);
}
generated quantities{
real s;
s=(p*(1-p))^.5;
}n=10
y=sample(0:1,n,replace=T)
summary(y)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 0.00 0.00 0.30 0.75 1.00par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')
data=list(N=n,y=y)
mdl=cmdstan_model('./ex1-1.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -7.53655
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 4 -6.10864 0.00573752 0.000129064 1 1 7
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -6.11
## p 0.30
## s 0.46mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -8.13 -7.86 0.68 0.31 -9.53 -7.64 1.00 916 957
## p 0.33 0.32 0.13 0.13 0.14 0.55 1.00 729 740
## s 0.45 0.47 0.05 0.04 0.34 0.50 1.00 719 740
geometric dist.
The event with probablity p occur after y trials(failure)
y~Ge(p), y[0,Inf), P(Y=y)=p(1−p)^y
E[y]=1/p, V[y]=(1-p)/p^2ex1-2.stan
data{
int N;
array[N] int y;
}
parameters{
real<lower=0,upper=1> p;
}
model{
for (n in 1:N) {
target+=log(p)+y[n]*log(1-p);
}
}
generated quantities{
real m;
m=1/p;
real s;
s=((1-p)/p^2)^.5;
}n=20
y=rgeom(n,0.3)
summary(y)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 0.75 1.00 2.20 4.00 8.00par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')
data=list(N=n,y=y)
mdl=cmdstan_model('./ex1-2.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -45.6597
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 4 -39.7495 0.00348 0.000336829 1 1 6
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -39.75
## p 0.31
## m 3.20
## s 2.65mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -41.80 -41.51 0.72 0.31 -43.22 -41.28 1.00 1117 1403
## p 0.32 0.32 0.06 0.06 0.23 0.42 1.00 759 1010
## m 3.22 3.13 0.61 0.55 2.36 4.39 1.00 759 1010
## s 2.67 2.58 0.62 0.56 1.80 3.85 1.00 759 1010
negative binomial dist.
The event with probablity p occur a times after y-1 trials
y~NB(a,p), y[a,Inf), P(Y=y)=p^a*(1−p)^(y-a-1)
E[y]=a*(1-p)/p, V[y]=a*(1-p)/p^2ex1-3.stan
data{
int N;
array[N] int y;
}
parameters{
real<lower=0,upper=1> p;
real<lower=1> a;
}
model{
y~neg_binomial(a,p);
}
generated quantities{
real m;
m=a*(1-p)/p;
real s;
s=(a*(1-p)/p^2)^.5;
int y1;
y1=neg_binomial_rng(a,p);
}n=50
a=3
y=rnbinom(n,a,0.3)
summary(y)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 4.00 6.50 7.38 9.00 23.00par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')
data=list(N=n,y=y)
mdl=cmdstan_model('./ex1-3.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -150.761
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 13 -142.64 0.00340739 0.00395503 0.9713 0.9713 16
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -142.64
## p 0.52
## a 3.81
## m 3.57
## s 2.63
## y1 7.00mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -144.10 -143.77 1.13 0.84 -146.39 -142.99 1.01 389 775
## p 0.63 0.61 0.17 0.20 0.37 0.93 1.01 248 453
## a 4.64 4.53 1.25 1.45 2.79 6.82 1.01 255 447
## m 2.75 2.82 1.37 1.54 0.50 4.89 1.01 264 455
## s 2.15 2.15 0.89 0.96 0.74 3.60 1.01 254 455
## y1 7.45 7.00 4.47 4.45 1.00 15.00 1.00 2076 1647
y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1') 
binomial dist.
The event with probablity p occur y times on n trials
y~B(n,p), y{0,1...n}, E[y]=np, V[y]=np(1-p)ex2-3.stan
data{
int N;
int n;
array[N] int y;
}
parameters{
real<lower=0,upper=0> p;
}
model{
y~binomial(n,p);
}
generated quantities{
real m;
m=n*p;
real s;
s=(n*p*(1-p))^.5;
int y1;
y1=binomial_rng(n,p);
}n=20
n0=10
y=rbinom(n,n0,0.3)
summary(y)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.00 2.00 3.00 2.95 4.00 5.00par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')
data=list(N=n,n=n0,y=y)
mdl=cmdstan_model('./ex2-3.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -136.737
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 4 -121.314 0.00270906 0.000861894 1 1 6
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -121.31
## p 0.30
## m 2.95
## s 1.44
## y1 2.00mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -123.34 -123.11 0.63 0.30 -124.53 -122.88 1.00 1130 1339
## p 0.30 0.30 0.03 0.03 0.25 0.35 1.00 879 1122
## m 2.98 2.97 0.31 0.32 2.49 3.49 1.00 879 1122
## s 1.44 1.45 0.04 0.04 1.37 1.51 1.00 879 1122
## y1 2.95 3.00 1.44 1.48 1.00 6.00 1.00 1763 1701
y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1') 
ex2-31.stan
trial n is varied from each sample
data{
int N;
array[N] int n;
array[N] int y;
}
parameters{
real<lower=0> p;
}
model{
y~binomial(n,p);
}n=20
n0=floor(runif(n,1,10))
y=rbinom(n,n0,0.3)
data=list(N=n,n=n0,y=y)
mdl=cmdstan_model('./ex2-31.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -49.9164
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 3 -49.3093 0.00401524 0.000319669 0.9516 0.9516 5
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -49.31
## p 0.32mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -51.35 -51.09 0.69 0.35 -52.78 -50.84 1.00 1144 1428
## p 0.32 0.32 0.05 0.05 0.24 0.41 1.00 809 1124
Poisson dist.
The event occur l times in unit range, the event occur y times.
y~Po(l), y{0,1...Inf}, E[y]=l, V[y]=lex1-4.stan
data{
int N;
array[N] int y;
}
parameters{
real<lower=0> l;
}
model{
y~poisson(l);
}
generated quantities{
real s;
s=l^.5;
int y1;
y1=poisson_rng(l);
}n=20
y=rpois(n,3)
summary(y)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.0 2.0 3.0 3.2 4.0 7.0par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')
data=list(N=n,y=y)
mdl=cmdstan_model('./ex1-4.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -22.8192
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 6 10.4417 0.000665407 0.000481145 1 1 8
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ 10.44
## l 3.20
## s 1.79
## y1 2.00mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ 11.15 11.40 0.63 0.29 9.84 11.61 1.01 1012 1316
## l 3.27 3.26 0.39 0.39 2.65 3.93 1.00 803 1224
## s 1.81 1.81 0.11 0.11 1.63 1.98 1.00 803 1224
## y1 3.28 3.00 1.86 1.48 1.00 7.00 1.00 1888 1944
y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1') 
exponential dist.
The event occur l times in unit range, the event take time y to occur
y~ex(l), y>0, E[y]=1/l, V[y]=1/l^2ex1-5.stan
data{
int N;
vector[N] y;
}
parameters{
real<lower=0> l;
}
model{
y~exponential(l);
}
generated quantities{
real m;
m=1/l;
real s;
s=1/l;
real y1;
y1=exponential_rng(l);
}n=20
y=rexp(n,2)
summary(y)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.007 0.095 0.378 0.448 0.537 1.920par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')
data=list(N=n,y=y)
mdl=cmdstan_model('./ex1-5.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -34.0293
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 5 -3.94689 0.000679157 1.92585e-05 1 1 9
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -3.95
## l 2.23
## m 0.45
## s 0.45
## y1 0.44mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -3.62 -3.34 0.71 0.30 -5.07 -3.12 1.00 846 1186
## l 2.35 2.30 0.51 0.49 1.57 3.27 1.00 741 1017
## m 0.45 0.44 0.10 0.09 0.31 0.64 1.00 741 1017
## s 0.45 0.44 0.10 0.09 0.31 0.64 1.00 741 1017
## y1 0.45 0.30 0.47 0.32 0.02 1.36 1.00 1853 1860
y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1') 
gamma dist.
The event occur l times in unit range, the event take time y to occur a times
y~Ga(a,l), y>0, E[y]=a/l, V[y]=a/l^2ex1-6.stan
data{
int N;
vector[N] y;
}
parameters{
real<lower=0> a;
real<lower=0> l;
}
model{
y~gamma(a,l);
}n=20
y=rgamma(n,2,1)
summary(y)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.41 0.79 1.60 2.03 2.55 6.28par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')
data=list(N=n,y=y)
mdl=cmdstan_model('./ex1-6.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -141.353
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 12 -32.3719 0.000120154 5.75695e-06 1 1 13
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -32.37
## a 1.83
## l 0.90
## m 2.03
## s 1.50
## y1 1.10mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -32.73 -32.45 0.98 0.78 -34.69 -31.75 1.01 669 1074
## a 2.12 2.07 0.59 0.59 1.23 3.10 1.01 617 728
## l 1.07 1.04 0.33 0.32 0.57 1.64 1.01 618 511
## m 2.03 2.00 0.33 0.32 1.56 2.61 1.00 2163 1619
## s 1.43 1.38 0.32 0.27 1.03 2.05 1.01 794 844
## y1 2.03 1.71 1.46 1.22 0.33 4.77 1.00 1849 1735
y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1') 
log normal dist.
logalithm of variable follows normal distribution
log y~N(m,s), y>0
E[y]=exp(m+s^2)
V[y]=exp(2*m+s^2)*(exp(s^2)-1)ex1-7.stan
data{
int N;
vector[N]<lower=0> y;
}
parameters{
real m0;
real<lower=0> s0;
}
model{
y~lognormal(m0,s0);
}n=30
y=rlnorm(n,0.5,1)
summary(y)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.2 1.2 1.6 3.2 3.1 35.6par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')
data=list(N=n,y=y)
mdl=cmdstan_model('./ex1-7.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -55.1291
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 8 -40.1814 0.00103138 0.00017249 1 1 11
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -40.18
## m0 0.56
## s0 0.92
## m 4.11
## s 3.12
## y1 2.01mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.3 seconds.seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -41.29 -41.01 1.02 0.77 -43.26 -40.31 1.00 987 1058
## m0 0.57 0.57 0.19 0.18 0.26 0.86 1.00 1282 1251
## s0 0.99 0.97 0.13 0.13 0.79 1.22 1.00 1894 1316
## m 5.05 4.54 1.99 1.38 2.97 8.70 1.00 1600 1334
## s 4.05 3.55 1.93 1.32 2.07 7.53 1.00 1668 1383
## y1 2.97 1.80 3.82 1.58 0.34 9.46 1.00 1945 2054
y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1') 
beta dist.
use as prior of binomial dist parameter p
event with probabilty p0 occur a-1 times, do not occur b-1 times in a+b-2 trials
p~Be(a,b), p[0,1], p=p0^(a-1)*(1-p0)^(b-1)
E[p]=a/(a+b), V[p]=ab/(a+b)^2/(a+b+1)ex1-8.stan
data{
int N;
vector[N] p;
}
parameters{
real<lower=0> a;
real<lower=0> b;
}
model{
p~beta(a,b);
}n=20
p=rbeta(n,1,2)
summary(p)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.006 0.225 0.293 0.301 0.339 0.607par(mfrow=c(1,2))
hist(p,main='p')
density(p) |> plot(main='p')
data=list(N=n,p=p)
mdl=cmdstan_model('./ex1-8.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -17.3024
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 7 9.06616 0.00223233 4.62945e-05 0.9944 0.9944 12
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ 9.07
## a 1.97
## b 4.71mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ 10.48 10.76 0.95 0.69 8.46 11.40 1.00 565 803
## a 2.23 2.16 0.61 0.60 1.32 3.31 1.00 652 748
## b 5.39 5.23 1.55 1.57 3.11 8.20 1.01 603 657
Dirichlet dist.
use as prior of categorical dist parameter p
p~dir(p0), p[0,1]ex1-9.stan
data {
int N;
int K;
matrix[N, K] p;
}
parameters {
simplex[K] p0;
}
model {
for(i in 1:N){
p[i]~dirichlet(p0);
}
}library(gtools)
n=20
p=rdirichlet(n,c(0.5,0.3,0.2))
summary(p)## V1 V2 V3
## Min. :0.003 Min. :0.000 Min. :0.000
## 1st Qu.:0.264 1st Qu.:0.096 1st Qu.:0.001
## Median :0.632 Median :0.314 Median :0.008
## Mean :0.530 Mean :0.331 Mean :0.139
## 3rd Qu.:0.762 3rd Qu.:0.483 3rd Qu.:0.117
## Max. :0.962 Max. :0.983 Max. :0.923boxplot(p)
data=list(N=n,K=ncol(p),p=p)
mdl=cmdstan_model('./ex1-9.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = 38.3768
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 5 56.0457 0.000115464 1.19815e-05 1 1 7
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ 56.05
## p0[1] 0.48
## p0[2] 0.34
## p0[3] 0.18mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ 51.50 51.84 1.07 0.70 49.37 52.47 1.00 1013 1180
## p0[1] 0.48 0.48 0.06 0.06 0.37 0.57 1.00 1866 1396
## p0[2] 0.34 0.34 0.06 0.06 0.25 0.43 1.00 1444 1147
## p0[3] 0.19 0.18 0.04 0.04 0.13 0.25 1.00 1847 1535
